近年来,随着我国一系列重点林业工程的实施,人工林面积已逐步接近全球人工林面积的30%(Xu et al.,2010)。发展人工林能够有效解决森林消耗、木材供应不足和森林生态环境被破坏等问题,而且在提高森林碳汇、维持区域生态平衡和加快森林恢复方面也具有重要作用(Payn et al.,2015;Szulecka et al.,2017)。当前,我国人工林经营已从单纯追求经济效益发展到注重多种效益相结合及可持续发展,其监测目的也从传统的森林资源利用转向保护生态环境、维持森林健康以及平衡森林碳储量等方面,高效和高精度的人工林森林资源监测越显重要(舒清态等,2005)。单木是构成森林的基本单元,其空间结构、生物物理和化学组分是森林资源调查、生态环境建模研究等所需的关键因子(李增元等,2016)。高分辨率精细遥感技术的发展,为人工林森林资源监测提供了有利条件,而快速、精确地获取单木信息可为森林资源监测与管理提供有效保障(Zhao et al.,2009;吴楠等,2017),合理利用各种单木分割算法,提高单木分割精度,可为精确获取单木空间位置、冠层结构以及掌握树木竞争和健康状况等提供重要的算法和技术支持。但传统的单木分割通常采用人工实测法,且仅能获得点上的数据,难以获取区域或更大尺度的数据(曹林等,2014);而利用高分辨率遥感手段提取单木信息可以提高森林调查效率,且能够在较高精度下保证单木信息的空间完整性和时间一致性(Wulder et al.,2010)。激光雷达(LiDAR, light detection and ranging)是一种快速获取地表目标三维空间信息的主动遥感技术,其通过发射激光脉冲到树木的枝叶等组分上,接收包含三维空间位置等信息的散射回波信号,并可生成点云数据,对精确估算冠层三维结构信息和单木分割具有一定的潜力和优势(刘鲁霞等,2014)。
目前,通过激光雷达数据进行单木分割的方法大致可分为2类:一是基于激光雷达点云数据生成栅格表面模型(CHM或DSM),在此基础上进行单木分割(Wang et al.,2004;Chen et al.,2006;Koch et al.,2006;Brandtberg, 2007;Popescu, 2007;刘清旺等,2008;2010;Yu et al.,2011),通常包括分水岭算法和多项式拟合法等;二是利用归一化后的LiDAR点云数据,根据LiDAR点云之间的空间结构关系和属性信息等进行聚类,直接进行单木分割(Alexander,2009;Reitberger et al.,2009;Li et al.,2012)。其中,分水岭算法、多项式拟合法和基于点云的距离判别聚类法应用较为广泛。Wang等(2004)在加拿大东北部采用分水岭算法对白云杉(Picea glauca)和少量冷杉(Abies fabri)、亚高山冷杉(Abies fargesii)进行单木分割,将自动探测的树冠图像与手动划定的树冠图像进行对比,结果发现75.6%的像元是相同的;Chen等(2006)在美国加利福利亚采用分水岭算法对阔叶林进行单木分割,并利用交叉验证方法进行精度评价,结果发现在阔叶林中单木分割的绝对准确率可达到64.1%。Popescu(2007)在美国东部采用多项式拟合法对松树(Pinus)林进行单木分割,结果表明多项式拟合法在该地区对上层和中上层单木探测效果较好,对下层单木的探测率不高(仅有28%);霍达等(2015)在内蒙古采用局部最大值法探测树冠,并用多项式拟合法拟合树冠边界,对以白桦(Betula platyphylla)为主的天然次生林进行单木分割,白桦冠幅的探测率达到80.5%,树高的平均误差为1.42 m,平均绝对误差为2.33 m,冠幅的平均误差为0.28 m,平均绝对误差为0.99 m。基于点云的距离判别聚类法的基本原理为利用栅格化的冠层高度模型(CHM)先探测树冠顶点,再用分水岭算法探测树冠边界。Alexander(2009)在芬兰南部采用基于点云的距离判别聚类法对松树、桦树(Betula)、云杉、花楸(Sorbus pohuashanensis)和毛白杨(Populus tomentosa)进行单木分割,结果发现75%的参考单木可以被识别;Reitberger等(2009)在德国东南部采用基于点云的距离判别聚类法对高山云杉林、山地混交林和云杉林进行单木分割,结果表明该三维分割方法能够在林分中低层检测小的单木,相对基于分水岭算法的单木分割方法精度最多可提高12%,单木分割准确率达到66%;Li等(2012)在美国加利福利亚内华达山脉采用基于点云的距离判别聚类法对混合针叶林进行单木分割,单木的探测率达到86%,分割准确率达到94%,整体精度F值达到0.90。
然而,以上对单木分割的研究大多只探讨了某一种单木分割方法的效果,并未在同一研究区对多种单木分割方法及其输入参数敏感性进行分析和比较;同时,以往的研究区大多位于亚寒带和温带森林,且单木分割对象多为天然林或次生林,基于激光雷达数据并针对亚热带人工林的单木分割相关报道并不多见。鉴于此,本研究以我国东南部亚热带人工林树种为对象,基于机载激光雷达点云数据,选取9个典型样地(3种不同林型,共378株),结合地面实测数据和目视解译方法,比较分水岭算法、四次多项式拟合法和基于点云的距离判别聚类法的单木分割精度,并对3种方法进行单木分割效果的敏感性分析,旨在:1)探索3种LiDAR单木分割方法对3种不同林型4个典型树种共378株森林单木进行分割,并通过计算单木探测率、准确率和F得分[F= 2(r×p)/(r+p)]评价分割精度;2)分析和比较不同林型内3种单木分割方法的单木分割精度;3)通过改变冠层高度模型(CHM)及调整基于点云的距离判别聚类法的距离阈值,分别对3种单木分割方法进行敏感性测试,探索不同算法参数对不同林型的敏感性差异。
1 材料与方法1.1 研究区概况及地面实测数据研究区位于江苏省常熟市虞山国营林场(120.70°E,31.67°N),所处区域为典型的北亚热带季风气候区,气候温和,年平均气温15.4 ℃,降水集中在6—9月,年平均降雨量约1 047.7 mm,年平均相对湿度80%左右。在研究区宝岩生态园内选择9块人工林样地(单木378株,样地大小30 m×30 m),主要树种包括水杉(Metasequoia glyptostroboides)、湿地松(Pinus elliottii)、香樟(Cinnamomum camphora)、杨梅(Myrica rubra)等。研究区树种及样地详细信息汇总如表 1、2所示。
表 1 研究区样地信息汇总Tab.1 Summary of samples in the study site表 1 研究区样地信息汇总
Tab.1 Summary of samples in the study site
样地编号No.密度Density/(tree·hm-2)类型Types平均树高Mean height/m平均胸径Mean DBH/cm样地1 Sample 1256水杉Metasequoia glyptostroboides16.8636.69样地2 Sample 2300混交林Mixed forests10.1422.96样地3 Sample 3344水杉Metasequoia glyptostroboides22.4144.53样地4 Sample 4478混交林Mixed forests10.0022.79样地5 Sample 5511混交林Mixed forests10.0023.22样地6 Sample 6544杨梅Myrica rubra7.3820.14样地7 Sample 7567混交林Mixed forests9.5521.00样地8 Sample 8633混交林Mixed forests7.8619.00样地9 Sample 9667杨梅Myrica rubra5.5412.03 表 2 研究区主要树种信息汇总Tab.2 Summary of forest characteristics for the four main tree species in the study site表 2 研究区主要树种信息汇总
Tab.2 Summary of forest characteristics for the four main tree species in the study site
林木参数Forest characteristics水杉M. glyptostroboides湿地松P. elliottii香樟C. camphora杨梅M. rubra平均值Mean标准差SD平均值Mean标准差SD平均值Mean标准差SD平均值Mean标准差SD树高Tree height/m18.924.2012.782.2412.704.326.201.87冠幅长轴Crown long axis/m5.841.546.761.618.182.677.221.89冠幅短轴Crown short axis/m5.291.514.951.315.552.095.991.67 1.2 遥感数据本研究采用Riegl LMS-Q680i全波形激光传感器获取LiDAR数据,数据采集时间为2013年8月17日。遥感平台飞行高度为900 m,旁向重叠度≥60%;传感器记录完整的激光脉冲返回波形信息,其时间采样间隔为1 ns,点云数据通过波形分解进行提取。激光发射脉冲所在波段为近红外(波长为1 550 nm),脉冲发射频率为360 kHz,最大扫描角为±30°;地面光斑直径约为45 cm,平均脉冲间距约为0.49 m,点云平均密度为8.74 pts·m-2。
1.3 研究路线地面实测数据包括样地中单木的位置、胸径、树高、冠幅、树种以及样地的定位信息,对区域LiDAR点云数据进行去噪,使用自然邻近法(natural neighbor)插值生成DSM和DEM,二者求差得到CHM,对CHM进行平滑处理后,利用分水岭算法和四次多项式拟合法进行单木分割;对于基于点云的距离判别聚类法,将离散的LiDAR点云数据借助DEM进行高度信息归一化,然后进行单木分割;最后,利用地面实测数据及部分目视解译数据分别对3种方法的单木分割结果进行精度验证和评价。机载激光雷达人工林单木分割方法比较和精度分析技术路线如图 1所示。
图 1 机载激光雷达人工林单木分割方法比较和精度分析技术路线Fig. 1 The technical route of comparisons and accuracy analysis of individual tree segmentation in plantation forests using airborne LiDAR data 1.4 分割方法1.4.1 分水岭算法分水岭算法是一种类似模拟浸水过程的数学形态学算法,将栅格图像CHM像元中的灰度值对应高度,然后将灰度值取反。图像中存在一些局部极小值点,假设在每个极小值点处打一个洞,从极小值处注水,水从图像中极小值点开始逐渐浸没,最终形成一个个积水盆。随着积水盆内的水位升高,不同极小值处区域的水面将汇聚到一起,此时在2个积水盆间建立一道堤坝;过程结束后,每个极小值点所在区域都被对应的堤坝包围,这些堤坝的集合便形成分水岭。堤坝是图像分割的边界,积水盆即为分割区域(Meyer et al., 1990;Popescu et al., 2004;Najman et al., 2005)。本研究首先通过分水岭算法确定树冠边界,然后在每个分割出来的边界中寻找局部最大值作为树顶(Walsworth et al., 1999;Wulder et al., 2000)。
1.4.2 四次多项式拟合法运用四次多项式拟合法的基本流程是先利用动态大小圆形滤波器(Popescu et al., 2004)在CHM中探测树冠顶点,基于每个树顶,根据单木和林分的形状结构进行四次多项式拟合(Popescu et al., 2003),确定树冠范围(图 2)。
图 2 四次多项式拟合法原理Fig. 2 Principle of polynomial fittingd1:东西向冠幅East-west canopy length;d2:南北方向冠幅North-south canopy length;d3:拟合出的冠幅Fitting canopy length;d4:实际冠幅Actual canopy length.通过冠层高度表面模型的栅格点拟合四次多项式,每个树顶进行2个方向的2次拟合,取平均值作为单木冠幅半径(霍达等,2015):
$f\left(x \right) = a{x^4} + b{x^3} + c{x^2} + dx + e;$ (1) $\mathit{f'}\left({{x_\mathit{n}}} \right) = 0{\rm{\;\;\;\;}}\mathit{n = }{\rm{1, 2;}}$ (2) $\mathit{f''}\left({{x_\mathit{n}}} \right) < 0{\rm{\;\;\;\;\;\;}}n = 1, 2.$ (3)四次多项式拟合的树冠曲线为式(1),式中x为拟合曲线的水平距离,a、b、c、d、e为四次多项式的拟合系数;式(2)中当一阶导数为0时,为拟合曲线的极值点;式(3)中当二阶导数小于0时,曲线为凹,为拟合曲线的极小值点,找到满足式(2)及式(3)的x1和x2,x1和x2为拟合窗口中凹点的位置,估测冠幅等于x1与x2的差值。
1.4.3 基于点云的距离判别聚类法基于点云的距离判别聚类法利用树木之间存在一定距离这一规律,假设归一化点云数据中最高点为树的顶点,然后从该点开始进行距离迭代判断,进行单木分割。其具体原理及过程如图 3所示。
图 3 基于点云的距离判别聚类方法单木分割过程Fig. 3 Process of point cloud-based cluster segmentationa:样地归一化后的点云A sample of normalized point cloud;b:基于点云的距离判别聚类法的判别原理Principle diagram of point cloud-based cluster segmentation;c:分割结果Segmentation result.点云按照从最高点到最低点的顺序进行分割,设置一个阈值,将离目标单木顶点水平距离大于阈值的点云排除出目标单木的点云集群,离目标单木顶点水平距离小于阈值的点云设置一个法则进行筛选提取。如图 3所示,目标树是树1,点A是树1的顶点。首先,点B是树2的顶点,排除出树1的点云集群,因为d1大于指定的阈值。接下来,d2小于指定的阈值,但点C属于树1,因为d2小于d3。点D不属于目标树1,属于树2,排除出目标树1的点云集群,因为d4大于d5。以此迭代判断,逐渐分割出单木(Li et al., 2012)。
1.5 精度验证指标采用3种不同方法提取的冠幅需要与地面实测和目视解译数据比较,并进行精度验证和评价。本研究采用以下3个指标进行衡量(Goutte et al., 2005;Sokolova et al., 2006):
$r = \frac{{{N_{\rm{t}}}}}{{{N_{\rm{t}}} + {N_{\rm{o}}}}};$ (4) $p = \frac{{{N_{\rm{t}}}}}{{{N_{\rm{t}}} + {N_{\rm{c}}}}};$ (5) $F = \frac{{2\left({\mathit{r} \times \mathit{p}} \right)}}{{r + p}}。$ (6)式中:r代表冠幅探测率,p代表探测出冠幅准确率,由r与p计算得出F得分;Nt、No、Nc分别代表正确分割、漏分和过度分割的数量。
2 结果与分析图 4为3种单木分割方法的分割效果,图 5为单木分割效果及其精度验证。
图 4 分割效果Fig. 4 Segmentation results 图 5 单木分割效果及其精度验证Fig. 5 Individual tree segmentation result and its accuracy assessment本研究总体F得分从高到低分别是基于点云的距离判别聚类法、分水岭算法和四次多项式拟合法,分别为0.83、0.81和0.76(表 3)。
表 3 激光雷达人工林单木分割F得分验证Tab.3 Accuracy assessments of individual trees segmentation based on airborne LiDAR data表 3 激光雷达人工林单木分割F得分验证
Tab.3 Accuracy assessments of individual trees segmentation based on airborne LiDAR data
精度验证指标Index for accuracy assessments分水岭算法The watershed algorithm多项式拟合法Polynomial fitting基于点云的距离判别聚类法Point cloud-based cluster segmentationr0.900.650.81p0.730.910.85F0.810.760.83图 6为分水岭算法和四次多项式拟合法的敏感性分析,分别选取0.3 m×0.3 m、0.5 m×0.5 m和0.7 m×0.7 m分辨率的CHM进行单木冠幅提取,并验证分割精度。
图 6 分水岭算法和多项式拟合法的敏感性分析Fig. 6 The sensitivity analysis of the watershed algorithm and polynomial fitting methods当CHM分辨率为0.3 m×0.3 m、0.5 m×0.5 m和0.7 m×0.7 m时,分水岭算法单木分割的F分别为0.74、0.84和0.60,四次多项式拟合法单木分割的F分别为0.69、0.75和0.45(表 4)。
表 4 分水岭算法与四次多项式拟合法单木分割精度敏感性分析统计Tab.4 The statistics of sensitivity analysis of the watershed algorithm and polynomial fitting表 4 分水岭算法与四次多项式拟合法单木分割精度敏感性分析统计
Tab.4 The statistics of sensitivity analysis of the watershed algorithm and polynomial fitting
分割方法Segmentation methods指标IndexCHM分辨率CHM resolutions/m0.3 m×0.3 m0.5 m×0.5 m0.7 m×0.7 m 分水岭算法The watershed algorithmr0.930.920.43n0.610.781.00F0.740.840.60四次多项式拟合法Polynomial fittingr0.790.630.29n0.620.941.00F0.690.750.45选取距离阈值,对基于点云的距离判别聚类法进行敏感性分析。对每块样地,阈值分别设为该样地最大冠幅半径、平均冠幅半径和最小冠幅半径。
由表 5可知,对于3种不同林型,当阈值设为样地平均冠幅半径时(阈值分别为2.93、2.51和2.23 m),分割精度最高(F分别为0.83、0.89和0.94)。
表 5 基于点云的距离判别聚类法单木分割精度敏感性分析统计Tab.5 The statistics of sensitivity analysis of the distance discriminant clustering method表 5 基于点云的距离判别聚类法单木分割精度敏感性分析统计
Tab.5 The statistics of sensitivity analysis of the distance discriminant clustering method
林型Forest types阈值(冠幅半径, m)Threshold value(crown radius, m)NtNoNcrpF简单林型Simple forest type最大Maximum=3.2027410.870.960.91平均Mean=2.2329220.940.940.94最小Minimum=1.7031051.000.860.92中等复杂林型Moderately complex forest type最大Maximum=5.1034920.790.940.86平均Mean=2.5137630.860.930.89最小Minimum=1.20394150.910.720.80复杂林型Complex forest type最大Maximum=5.70381960.670.860.75平均Mean=2.93471080.820.850.83最小Minimum=1.00507210.880.700.78 3 讨论本研究选取CHM分辨率作为分水岭算法和四次多项式拟合法的关键参数进行敏感性分析, 选取距离阈值作为基于点云的距离判别聚类法的关键参数进行敏感性分析, 结果发现:1)当CHM分辨率为0.5 m×0.5 m时, 分水岭算法和四次多项式拟合法的分割精度最高。CHM分辨率过小, 分割时会导致多分现象增多; CHM分辨率过大, 分割时会导致漏分现象增多。多分和漏分都会影响分割精度, 这是因为当分辨率过低或过高时, 插值生成的CHM与真实的林分冠层高度契合度不高, 分辨率过低会出现许多不符合真实林分冠层高度的凸和凹, 产生过度分割; 分辨率过高则过于平滑, 不能详细体现真实林分冠层高度, 产生漏分(Sokolova et al., 2006)。2)当阈值近似平均冠幅半径时, 基于点云的距离判别聚类法分割精度最高。所选距离阈值小于该林分平均冠幅半径, 过度分割现象增多, r值偏低; 所选距离阈值大于该林分平均冠幅半径, 漏分现象增多, p值偏低。
对运用分水岭算法进行单木分割的结果进行精度评价后发现, 其p值较低, 表明该方法过度分割的数量较多, 过度分割是影响该方法分割精度的主要因素; 对运用多项式拟合法进行单木分割的结果进行精度评价后发现, 其r值较低, 表明该方法的漏分较多, 漏分是影响该方法分割精度的主要因素。Larsen等(2011)使用欧洲3个国家6个不同位置、不同森林类型的数据, 对6种单木树冠提取方法进行评估, 包括局部最大值法(local maxima detection)、谷地追踪法(valley following)、模板匹配(template matching)等, 但其并未对各分割方法选取关键参数进行敏感性分析; Kaartinen等(2012)将多种单木提取方法汇总到一起进行综合对比, 包括基于栅格CHM的多种单木提取方法和基于点云的聚类分割方法, 结果发现所有单木提取方法使用效果均较好, 对单木树高的提取精度均优于0.5 m, 然而此种方法的研究区并不统一。今后在利用分水岭算法进行单木分割研究时, 可以考虑运用标记控制分水岭算法, 其适应性更强, 可有效减少单木分割时的过度分割现象, 如Culvenor (2002)提出的树木识别和描绘算法TIDA (tree identification and delineation algorithm)等。
本研究发现,对于不同类型林分,需要调整各单木分割算法的关键参数使其更适应此林分,不同算法面对不同类型林分各有优劣。基于LiDAR的单木分割算法大多基于单木的结构特征,本研究获取的LiDAR点云平均密度为8.74 pts·m-2,若提高点云密度(如通过低空多旋翼无人机获取极高密度点云),应能使其更好反映单木和林分的结构特征,有利于各算法单木分割精度的提高;此外,若结合LiDAR点云数据与光谱特征,也将具有提高单木分割精度的潜力(申鑫等,2016)。本研究比较的分水岭算法、四次多项式拟合法和基于点云的距离判别聚类法分别是基于栅格数据与点云的单木分割算法,近期也有研究结合栅格数据与点云数据进行的单木分割算法(Dalponte et al., 2016),今后可尝试比较